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ABSTRACT: 



The response of a plane Poiseuille flow to disturbances of 
various initial wavenumbers and amplitudes is investigated by 
numerically integrating the equation of motion. It is shown that 
for ver y lQW _amul itude disturbances the numerical integration scheme 
yields results that are consistent with those predictable from_ 
linear theory. It is also shown that because of non-linear interactions 
a growing unstable disturbance excites higher wavenumber modes which 
have the same frequency, or phase velocity, as the primary mode. 

For very low amplitude disturbances these spontaneously generated 
higher wavenumber modes have a strong resemblance to certain modes 
computed from the linear Orr-Sommerfeld equation. 

In general it is found that the disturbance is dominated for a 
long time by the primary mode and that there is little alteration of 
the original parabolic mean velocity profile. There is evidence of 
the existence of an energy equilibrium state which is common to all 
finite -amplitude disturbances despite their initial wavenumbers. 

This equilibrium energy level is roughly 3-5$ of the energy in the 
mean flow which is an order of magnitude higher than the equilibrium 
value predicted by existing non-linear theories. 
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1. Introduction 



In this paper we present some results from a numerical investigation 
of the phenomena of stability, transition and turbulence in incompressible 
channel flows. Ideally, the goal of such an investigation would be to 
numerically integrate the Navier-Stokes equations and thereby be able 
to follow the growth of some unstable disturbance which perturbs an 
initial laminar Poiseuille flow. It would then be possible to observe 
in detail how this growing disturbance eventually produces a fully 
developed, stationary turbulent flow. Presumably, if one were to repeat 
such a calcxolation for a large number of initial disturbances then the 
final turbulent flow could be described by an ensemble average over all 
the initial conditions . 

The results of such a computation would represent an exact (to 
within the limits of the numerical model) solution to the turbulence 
problem. Since turbulence is an inherently three-dimensional phenomenon, 
these computations would have to be carried out on a three-dimensional 
grid; unfortunately, the computational grid would have to be large enough 
to contain the large "eddies” representing the initial disturbance and 
still be able to resolve the smallest energy dissipation length scales 
of the final turbulence . Even at low Reynolds numbers the ratio of the 
larger length scale to the smallest scale is several orders of magnitude. 
The number of grid points is proportional to the cube of this ratio; 
hence, a complete three-dimensional solution is impractical, if not 
impossible, despite the large storage capacity of modern computers. 

(Emmons 1970) To overcome this impractically large storage requirement 
we have drastically simplified the problem by treating the flow on a 
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two-dimensional basis. It is recognized that a two-dimensional treatment 
can not adequately represent true physical turbulence. Never-the-less 
the two-dimensional problem can be handled without further approximation 
and, as will be pointed out in the next section, two-dimensional 
solutions are of interest in their own right. 

2 . Background 

During the past decade considerable effort was devoted to developing 
theories for the response of plane Poiseuille flows to finite -amplitude 
disturbances. The most significant contribution of these theories 
(in particular, those of Stuart and Watson (i960) and Reynolds and 
Potter (1967)) was in showing that the square of the amplitude ( |A^ | ) 
of an initially infinitesimal unstable disturbance is governed by an 
equation of the form 

d I A 1 I 2 p u 

it = *1 l A l ! +k 2 1*1 1 (1) 

Equation 1, which can be derived from the Navier-Stokes equations, is 
an approximation valid in a region in the ( or-R e ) -plane which is close to 
the "neutral curve". Here R g denotes the Reynolds number and a the 
wave number of the disturbance. For disturbances with very small 
amplitudes the second term on the right of equation 1 must be neglible; 
hence, the square of the disturbance amplitude has an exponential 
growth and the constant, k^, is related to the exponential amplification 
factor of the linear theory (i.e., k^ = -2^1 |3_^). Of particular interest 
is the case of k^ > 0 for which the flow is unstable to small disturbances. 
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If k^ <0 then it is possible that the higher order terms will 
eventually balance the leading term and the amplitude growth is limited 
such that 



a M 2 

ST — - 0 

dt 

K I 2 - 



(ib) 



as 



t -♦ + oo 



Under such conditions a "supercritical equilibrium state" is 
said to exist. On the other hand, if > 0 and k^ > 0 then no 
positive limiting value for |A^ | is possible according to equation 1 . 

It is possible to evaluate k^ based on Stuart’s (i960) theory and 
the necessary calculations have been carried out by Pekeris and 
Shkoller (1967) and Reynolds and Potter (1967). These authors have 
found that under certain conditions k^ < 0 for some unstable disturbances. 
Both Pekeris and Shkoller and Reynolds and Potter have found that for 
regions near the lower branch of the "neutral curve" k^ < 0 while 
k^ > 0 near the upper branch. Hence, one can conclude that disturbances 
whose wavelength corresponds to a point near the lower branch can grow 
according to the non-linear theory and reach a supercritical equilibrium 
state. On the other hand, for unstable disturbances whose wavelength 
corresponds to a point in the (a, R g ) plane near the upper branch the 
non-linear theory does not seem to be applicable. At least the existing 
theories give us no information about the growth of such a disturbance 
beyond the linear range of amplitudes. 
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Such "supercritical equilibrium states" are not observed experimentally 
in channel flows; instead the initial instabilities lead immediately 
to three dimensional turbulent flows. It seems, that no experimental 
verification of this aspect of finite -amplitude theory is possible. 

A computer simulation of a two dimensional channel flow is, perhaps, 
the best way of investigating the two dimensional response of a plane 
Poiseuille flow to an initial unstable disturbance. A computer 
simulation offers, in addition to its two dimensionality, the advantage 
of providing an "exact" solution to the equations of motion ("exact" 
to within the limits of a finite difference representation) without 
introducing any approximations concerning the relative magnitudes of 
the various modes of which the disturbance is composed. 

The present paper presents the results of a numerical investigation 
of the response of a plane Poiseuille flow to various initial disturbances. 
An attempt has been made to present the results in a manner that will 
facilitate comparison with the existing non-linear theory. The 
difficulty in making such comparisons is due to the fact that the 
theories mentioned concern the growth of an unstable disturbance with 
z ercu ini t i al amplitude (at time t = - « ) which grows at first according 
to linear theory and finally reaches an amplitude of such magnitude 
that the non-linear equations are appropriate. In a numerical simulation, 
with limited computation time available, the initial disturbance must 
have some finite amplitude. The growth rates predicted by the finite- 
amplitude theories are very small and the time required to follow 
the growth of a very small disturbance to the non-linear range consumes 
enormous amounts of computer time. For this reason the numerical 
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investigator is tempted to start his calculations by using a disturbance 
whose amplitude ,is.,pf_such~mgnitui^jbh_ at the e ntire linear range can 
be b ypassed. It will be pointed out in Section 6 that there seem to be 
significant qualitative differences in the flow when large initial 
amplitudes are used as compared with those for which very small initial 
amplitudes are used in the calculations. 

Another aspect of the present investigation was to determine under 
what circumstances, if any, a two dimensional representation of the 
flow could be made to resemble in some respects three dimensional 
turbulent flows . 

It is clear that the "supercritical equilibrium states" in which 
the unstable disturbance mode remains dominant do not resemble 
turbulence in which the disturbance energy is distributed over a 
wide spectral range. AJ_so in a "supercritical equilibrium state" 
the resultant mean velocity profile differs only slightly from the 
parabolic laminar profile while for turbulent flows the mean velocity 
profiles bear little resemblance to the original laminar one. It 
is recognized that there are significant differences between two 
dimensional and true three dimensional turbulence. This is particularly 
evident in the energy spectrum and mechanism for the transfer of 
energy and vorticity between the various spectral components (Lilly (1968) 
and Kraichnan (1967 ) ) • 

However, could it be possible that for some choice of initial 
disturbance mode shapes and amplitudes the structure of the resultant 
two dimensional flow could lie somewhere between the (relative) 
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simplicity of a supercritical equilibrium state and the "complete 
chaos" characterizing a true turbulent flow? This aspect of the 
present paper will be discussed in Section 6. 



3. Basic Equations 

Consider the two dimensional flow of an incompressible fluid 
between two parallel planes . The governing equations of motion are 



K _ *1 K +*! K = 2 r 

3t By 3x Bx By R e ^ 

and 

2 

C = v ¥ 



(2a) 



(2b) 



where £ and Y are the vorticity and stream function respectively with 
£ and Y defined by 



u 



BY . _ BY 

By ’ v 3x 



(3a) 



r _ Bv 3u 
c ~ Bx " By * 



(3b) 



The variables are assumed to be normalized on suitable characteristic 
lengths and velocities which are shown in Figure 1 along with other 
parameters of interest* The velocity components can be expressed as 



U = U(y) + u' ; V = V* ; £ = - ^ + C 

Y = 3/2 (| y 3 - y) + Y' 



(4) 



where the primed quantities which depend on x, y, t represent the 
instantaneous departure of the flow from the original laminar flow . 
Substituting equations h into equations 2 results in 
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(5a) 



C = V 2 Y 



(5b) 



The conditions of zero slip at the wall and constant mass flow rate 
imply that the stream function, Y ! , satisfy the boundary conditions 



By considering periodic solutions to equation 2 the boundary conditions 
in the streamwise direction are fixed by 



where the basic period is 2L. Equations 5& and 5b with the boundary 
conditions (6) and (7) are to be solved for Y T (x,y,t) for a given 
initial value of Y T (x,y,0). 

At this point it is convenient to define some parameters that are 
useful in describing the flow structure. The assumed periodicity 
implies homogeneity in the streamwise direction; hence, all averaging 
is done with respect to the streamwise coordinate, x. The average 
value of any quantity, q(x,y,t), is denoted by an overbar and obtained 
by integrating over one period of the basic wavelength. 



r = SY'/Sy =0 at y = ± 1. 



(6) 



y (x»y»t) = y ( x ± 2nL, y, t) 



(7a) 



C (x,y,t) = C' (x ± 2nL, y, t) 



(7b) 




( 8 ) 
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It should again be emphasized that in our representation the flow' is 
arbitrarily divided into a laminar and a disturbance part. Hence, 
the total streaiiMise velocity component is given by 

u(x,y,t) = U(y) +u’(x,y,t) (9) 

and the mean velocity is 

u(y,t) = U(y) + u r (y,t) (10) 



The difference between the mean and local values of any quantity is 
referred to as the turbulent part of that quantity and is denoted by 
a circumflex. Thus we have 



u (x,y,t) = u(x,y,t) - 


u(y,t) = u’ - u’ 


(Ha) 


v (x,y,t) = v’ (x,y,t) 




(lib) 


T (x,y,t) = V (x,y,t) 


- r" (y,t) 


(He) 


C (x,y,t) = C’ (x,y,t) 


- V (y,t) 


(lid) 



and 



A 

U - 



3Y 
3 y 



u’ 



A 

V 



ar 

3y 



A 

3x 



v* 



3T 

Sx 



= 0 



The turbulent kinetic energy is 



E (x,y,t) = \ [(- §) 




(lie) 

(Ilf) 

(Ug) 

(Hh) 

(12a) 
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The mean value is 



— 1 * * 

£ (y 9 t) = PL J E ( x jyj t ) dx 



(12b) 



and the normalized sum of turbulent kinetic energy over a domain extending 
from the upper to the lower wall and for one period of the basic 
disturbance is 



E(t) = 



m 1 * r k t 

-1 L -L 



L A 

E (x,y,t) dx! 



(12c) 



Similarly, the kinetic energy in the mean flow is 



E (y,t) = i [u(y) - 



(13a) 



and the normalized sum over the domain is 



E (t) = h J E (y,t) dy 

-1 



(13b) 



The normalized energy in the initial laminar flow is 



E^i | \ U 2 dy = .600 . 


(13c) 


-1 




can represent the periodic disturbance by means 


of a Fourier 


as follows : 




V (x,y,t) = E f n (y,t) e -1 nax 

n=-oo 


(14a) 


C (x,y,t) = s g n (y,t) e‘ inQK 

n=-co 


(14b) 
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with 



max 



f (y,t) = 57 J V (x,y,t) e 1 dx 
n 2L -L 



(15a) 



f = f * 
-n n 



g n (y,t) = ^ J V (x,y,t)e inQK dx ( 15 b) 

— L 

«-» - C 

where a - tt/L ? and superscript * denotes a complex conjugate. From 5b 

g = f " -(n a) 2 f (15c) 

n n n 

Notice that the mean value of the disturbance streamfunction and vorticity 
is simply 



r (y,t) = f Q (y,t) (l6) 

V (y,t) = g Q (y,t) (17) 

The mean value of the turbulent kinetic energy is given by 

S (y,t) = S E n (y,t) - E Q (y,t) (l8) 

n=-° o 

where each of the components of the energy spectrum function, E (y,t), 
is given by 



E n (y,t) 




(19) 
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The normalized total turbulent kinetic energy over the domain is 



XV a* 



E(t) = E E n (t) - E o (t) 



(19a) 



where 



1 



E n (t) =| J E n (y,t)dy 



(19b) 



4. The Numerical Model 

A rectangular grid of dimensions 2 x 2L is used with M x N grid 
points located as shown in Figure 2. For the calculations presented 
herein 64 x 201 grid points were used. The indices k, j , and i denote 
stations along x, y, and t respectively and 6x, 6y, 6t are the 
corresponding intervals . The value of the disturbance streamfunction 
r (x,y,t) is then denoted by Y™. T . The use of the primes will be 

IV j J 

discarded from this point unless necessary for clarity. Equations 
5a and 5b are expressed in finite difference form as follows 




3 



K+l.J K-l.J 



^K+l , J ~ ^K-l, J 
2 6x 



2 6t 



26x 



(20a) 





Y 1 ) 

t K-1,J / 



(20b) 
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2 

where A represents the modified form of the Laplacian differential 
operator of DuFort and Frankel (1958) defined as 



■ l - 1 r r l (r l+ 1 + + r 1 1 

K,J 6x 2 L 4+1,1 ' V4,J 4,JV 4-1, jJ 



+ 




°K, J+l 





+ c 



x 1 

K, J-l J 



( 21 ) 



and j is the total energy and mean square vorticity conservation 
form of the non-linear advection terms introduced by Arakawa ( 1966 ). 
i.e.. 



J K,J 3 { A y l K,J A x "K,J ■ A x f K, J A y "K, j} 

+ 5 K (ta =c c) 1 j (23) 

+ I { A x<£ V>K,J - A y«= ■ 



where A and A denote ordinary first central differences . 
x y 

A+l 

Equation 20a is solved explicity for £ at each time step from 

IV J J 

a a & -l 

the known values of T , A T > and T • The corresponding values 

iVjJ iv 5 J 

^_j_l 

of the Y are obtained from equation 20b. Equation 20b represents 
K 5 J 

a set of simultaneous algebraic equations to be solved for the Y 

K j J 

K = 1,2,... M; J = 1,2,... N at each time step. Because of the large 

number of mesh points used in this investigation, conventional 

techniques are not suitable and an alternative method is used. We 

i & 

express the variables T and T in terms of their discrete Fourier 

Iv 2 J I\ j J 

transforms as follows : 
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(24a) 



Y* = Z f 1 e _i 0 n,K 
KjJ 21=1 



M 

r SL i -i 0 n k 

Ek,j * Z g n,J e n,K 
’ n=l ’ 



(24b) 



with 



M 



f \ = i Z Y* e 1 0 n,K 
n,J M k=l K ’ J 



(24c) 



M 



^ = - Z r 1 e 1 0n > K 

T M ^ W T e 



3 n, J M 



k=l 



K, J 



(24d) 



where 



Q = 2 tt (n-l) (k-l) 
" M 



n ,K 



k,n = 1,2, . . .M 



Substituting equations 24 into equation 20b yields the following 



equation for the Fourier components of Y„ T : 

K, J 

* 2 
f TJ , - a f , + f T , = 6 g 
n,J+l n n,J n,J-l y n,J 



(25) 



n ~ 1>2, . . .M 
J = 1,2, ...N 



where 



a* = 2^1 - (6y/6x) 2 (cos 9^ - l)] 



The transforms and inverse transforms defined by equations 24 are 

computed by the fast Fourier algorithm of Cooley and Tukey (1965) while 

equations 25 are tri-diagonal and are easily solved for the f 1 s 

n $ J 

by Gauss elimination. 
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The assumed periodicity in x is automatically satisfied through 

equations 24. The boundary conditions at the walls (equation 6) are 

0 

solved by setting Y„ = Y„ , T = 0 for all times and the no slip 

Ji j X IV y iM 

jj 

condition is satisfied in the following way. We express T in 

j 

terms of Y^ , through a Taylor series expansion about a point on the 
K-ji 

wall. 



’X - + *1 Sy + A 2 6y2 + A 3 Sy3 + ••• 

Yr ^3 = YjA + 2A 1 6y + 4A 2 6y 2 + 8A 3 6y 3 + ... 
Y K^4 = 2^ Sy + 9A 2 6y 2 + 27A 3 6y 3 + 



(26) 



etc . 

The no slip condition (dY'/^y = 0) is satisfied if we require that 

£ 

A n = 0 which implies that (since Y^ , = 0) 

1 IS.,1 



and 



for all k 
order) at 



\Lf ^ - 2. yu ^ 

*K,2 2 *K,3 9 K,4 


(27a) 


w ^ = I f ^ _ i ^ 

t K,W-1 2 t K,N-2 9 K,N-3 


(27b) 


and 1. The Poisson equation (5b) is satisfied (to 
the walls by 


second 


£ K,1 = ^2 = ~^2 { 3 / 2 Y K,3 " ^ 9 Y K,4} 


(28a) 


£ K,N = ~^2 { 3//2 Y K,W-2 " Y K,N-3} 


(28b) 
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* 



The length of the time increment, 6t, is fixed by the semi -empirical 
stability limit 



6t = f 6x 6y/£(U + u) q 6y + v q ' 6 xJ 



( 29 ) 



where the velocity components (U+u') to v q ' are the maximum absolute 
magnitudes of velocity selected from a large sample of grid points. 
Numerical experimentation has shown that for values of f ^ .6 consistent 
results are obtained. Values of f near unity may lead to catastrophic 
instability in this explicit technique. 

Since equation 20a is not "self starting" (it requires three time 
levels), the computation procedure starts with the forward difference 
representation 



/ +1 . c 1 

r T 



Y t ^ 1 t _ Y S/+1 



i+i 1+1 

K, J ~ b K,J _ I r * K+1 , J ~ *K-1,J °K+1 , J " °K-1, J 

6t 2 r 2 6x " U J 26x 



+ / +1 + -L A 2 c^ +1 ) 
d K,J R ®K, Jj 



Y ^ _ w ^ r ^ r ^ 

i | 3 K+1,J _ t K- 1,J _ y °K+1 , J t - 



K-l, J 



26x 



2 6x 



(30) 



T A 1 a 2 

J K,J + R e 



A 2 c n 

wC, JJ 



Equation 30 is also used periodically during the calculations in order 



to suppress any instabilities associated with the central time differencing 



used in equation 20a. 



The computation procedure can now be outlined step-by-step. 

1. An initial distribution of ¥„° T is chosen and the corresponding 

K • J 

O * 

Q are computed from equation 20b. 

J 
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2. The time increment, 6t, is computed from equation 29. 

3. The values of are computed for all the interior points 

(j = 3,4,...N-2) by iteratively solving equation 30. 

4. The Fourier coefficients, are computed from equation 

24d, equation 25 is solved for the f 1 T and the corresponding 

q n, J 

values for the j found from equation 24a for the interior 
points (J = 3,4, . . .N-2) . 

5. The remaining grid points for J = 1,2,N-1,N) 

are computed from equations 27a, 27b, 28a ,b and 20b. 

6. For each time step the procedure in steps 3, 4, and 5 is 
repeated except the explicit central difference equation 20a 
is used in step 3 rather than equation 30 . 

7. Periodically, say every 50 time steps, a new time increment, 

5t, is computed (step 2) and a single time step using the 
implicit forward difference equation 30 is used. 

For all of the computations presented in this paper a grid with 
dimensions M = 64, N = 201 was used which required a computation time 
of approximately 17 seconds per time step. 

5. Linear Calculations 

In an investigation such as this it is desirable to make some 
comparison between the results of our numerical computations and some 
known solution to the Navier Stokes equations. The imposed restraint 
of two dimensionality precludes any comparison with experimental data 
and, of course, there are no exact solutions to the non-linear equations 
of motion. The only alternative is to consider the limiting case of 
infinitesimal disturbances for which the results from linear theory 
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are available. For this reason, a rather thorough investigation of 
the linear solutions to equations 5a and 5h was undertaken. For small 
amplitude disturbances the non-linear terms in the vorticity transport 
equation are neglected and the equations of motion become 



St 



+ u 



K' 

Sx 





(31a) 



C = V 2 Y* 



(31b) 



The general solution to equations 31 can be expressed in the form 



03 00 



r(x,y,t) = E E $ nm (y) e 

n=-co m=l 



-i(n ax - p t) 



nm 



(32a) 



co oo -i(n ax - p t) 

C (x,y,t) = E E e nm (y) e 

n=-oo m=l 



(32b) 



with 



e = $ - (n a) $ 

nm nm ' 7 nm 



(32c) 



a = 2 tt/x 



where X is the wavelength of the disturbance. 

Substituting equation 32 into 31 results in the well known Orr- 
Sommerfeld equation for $ nm (y)j i-e., 



$ 

nm 



2 

2(na) 

- (m > 2 o 



(na)^ 

T1 

+ U 



$ + i na R Tf u - — ' )(* 

nm e LA na A i 

$ 1 = 0 
nmj 



(33a) 
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with boundary conditions 



$> = $ =0 at y = - 1 . 

nm nm 



(33b) 



Solutions of 33a for $ nm (y) for given a, n, and R g can be made to satisfy 
the homogeneous boundary conditions (33b) only for the discrete eigen- 
values, p , 1 £ m £ ». For the given a, n, and R g we are usually 
interested in determining if there exist corresponding p^'s which 
have negative imaginary parts. The corresponding disturbance, ^^(y)? 
then grows exponentially with time and the flow is unstable to disturbances 
at the specified a , n, and R g . 

For the present analysis we are interested in solutions to equation 
20a; hence, we must consider the finite difference representation of 
the linearized equation. For small disturbances equation 20a becomes 

J l ...A . SL 



£+1 _ -- G -1 - Y ~ - C 

fa K,J b K,J K+1,J t K- 1,J b K+l,J b K-l,J 

26t 5 2 6x U J 2 6x 



+ VR e A 2 (J 



Again we will consider solutions of the form (32); hence, we have 



( 3*0 



Y ^ = ( § ) e 

*K,J v Vnm 



-i(na(k-l)6x - p i 6t) 



nm 



^K, J = ^®J^nm e 



-i (ncv(k-l) 6x - p SL 6t) 



( 35 ) 



nm 



Substituting into 3**- and 20b results in 



OH™ = W [W +3<i j) ] + ryy [ (l W nm + (0 j- i) 1 

nm nm R^Sy nm nm nm 



nm nm" 

( 36 a) 
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+ 0 ) 



(36b) 



<v 

nm 




<*j) 



nm 



nm 



where 



u> =2^(6y/6x) 2 cos na 6x^ 

2 

u) = -2^1 + (j~) (1 - cos na 6x)J 



a' = sin na 6x/n6x 



and 



p' = 2[l + l/R 0 6x 2 l cos p 6t + i 
nm L e J nm 



sin P 6t 
nm 

6t 



= V-1 



(37a) 

(37b) 

(37c) 



(37d) 



When expression 36b for (0 T ) is substituted into equation 36 a there 

d nm 

results a homogeneous set of linear algebraic equations for the (§ ) 

J nm 

For given values of a, n, R g , 6x, 6y, and 6t we wish to determine the 

eigenvalues and eigenvectors of the matrix of coefficients of the ($ ) 

J nm 

In particular, we are interested in those eigenvalues p for which the 

nm 

imaginary part of is negative and in the associated eigenvectors 

($ ) . The effects of 6x 3 6 y 3 and 6t on the linear solution are especially 

J nm 

important for the present application. 

The eigenvalues and eigenvectors have been computed for a number of 
Reynolds numbers and for a wide range of a's using the QR algorithm 
(Wilkinson 1965 or Fair 1971)* For more details on these calculations 
the reader is referred to 0 f Brien (1970) or to the report by Gawain and 
Clark (1971). In the latter report some of the eigenvalues (p^) and 
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eigenfunctions (^^(y)) are tabulated for the case a = 1.0, R g = 6667 . 

In the limit as 6x, 6t -* 0, the problem is the same as that solved by 
Thomas (1953) although the methods used are quite different from those 
employed by Thomas. For the case a = 1.0, R g = 6667 , and 6y = .01 (in 
Thomas' notation this corresponds to a = 1.0, R g = 10,000, 6y = .01) we 
have solved equations 36 for n = 1, 2, 3 5 and 4. Our solutions show 
that there is only one unstable eigenvalue and this unstable solution 
occurs for n = 1. We shall order the eigenvalues, (3^, so that for a 
given n, the m's are arranged according to the stability of the computed 
eigenvalues. Hence, for the case a = 1.0, R e = 6667 , is the unstable 
eigenvalue with the corresponding eigenfunction $qq(y)* Thomas gives a 
value of (3^ = . 3&53 - i *0055 and the corresponding mode shape, ^ n (y) 
is tabulated in his Table V. (Using our non-dimensionalization, is 
exactly 1 . 5 times the value given by Thomas . ) Our calculations yield 
the value (3^ = .3593 - i .004l and the corresponding mode shape, §^(y), 
is shown in Figure 3 for comparison with that of Thomas . 

Now consider the influence of the finite difference parameters 
6x, 6y, and 6t on the solutions to 36 . The effect of 6y on the growth 
rate (-^ (3^) i s illustrated in Figure 4. As 6y increases the flow 
becomes more and more stable and when 6y = . 05 , the phenomenon of linear 
instability has been completely obscured by the crude finite difference 
representation of the equations of motion. Indeed, calculations with 
iy = .05 for a wide range of Reynolds numbers and o's have shown no 
indication of linear instability. 

The effect of varying 6x was investigated at R g = 6667 , o = 1.0, 

6y = .01 and 6t < .02. Figure 5 shows the variation with 6x of the growth 
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rate (-J P 11 ) and phase velocity (^ 1X ) for the unstable mode . For 

6x < .8 the changes in the disturbance profile (§^(y)) where insignificant 

Perhaps the most surprising discovery in this investigation of the 
linear behavior of the discrete representation of the equations of motion 
was the significance of the time interval, 6t, on the disturbance growth 
rate (4 (3^). The effects of 6t on this parameter are shown in Figure 6. 
There are significant departures from the "exact 11 value of>J p for 
6t > .04 while the phase velocity and mode shape are not significantly 
changed. Hence, in addition to the usual stability constraint (29) on 
the time interval , 6t, it is necessary that 6t < .04 in the numerical 
integration scheme described in section 4. 

We can now summarize the results presented in the preceeding discussion 
Increasing the y increment , 6y, tends to make the flow more stable and 
for coase grids ( 6y > .05) the phenomenon of linear instability is 
completely obscured. The linear solution is comparatively insensitive 
to changes in 6x, at least, until the x-increment is approximately the 
same length as the channel half -width. The growth rate of the disturbance 
is strongly influenced by the 6t increment (due to the DuFort-Frankel 
representation of the diffusion term) and for large 6t the growth rate 
is much higher than that predicted by the "exact" linear solution. 

6. Numerical Integration of the Vorticity Equation 

(a) Small Amplitude Disturbances 

We now return to the numerical integration of the vorticity 
equation using the techniques described in Section 4. As a check on 
the validity of the solution we first consider the limiting case of 
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small amplitude disturbances and compare the results obtained with the 
known results from the linear solution as described in the preceeding 
section. To be more explicit, for a given Reynolds number and initial 
disturbance streamfunction, Y f (x,y, 0 ), are the results from the numerical 
integration consistent with the values of and ^^(y) predicted by 
linear theory? 

A series of calculations were made at R = 6667 and the initial 

e 

disturbance given by 

Y* (x,y, 0 ) = e jf^Cy-jO) e“ 1QK + f^(y,0) e +1 ax j- ( 38 ) 

where e is an amplitude factor. For e sufficiently small the non-linear 

terms in equation 5 a should be neglible and the numerical integration 

should give results which are in agreement with the linear theory. As 

a first test we choose a = 1.0, R = 6667 % € = .005 and the initial mode 

e 

shape identical to the unstable mode computed from linear theory; i.e., 
f-L (y»0) = $ 2.1 ( y )> where ^(y) was shown in Figure 3. The growth 

A 

rate of the disturbance kinetic energy, E (t), as defined by equation 12c 
is shown in Figure 7. The computed growth rate gives a value = -0050 

which is in good agreement with the value given by Thomas. From equation l4 
we can determine the phase, cp n (y,t), of any of the spectral components 
of 7 (x,y,t) by 



%(y»t) = tan ” 1 (4 ^n/j^f ) (39) 

The phase velocity, c , or rate of propagation, of any of the spectral 
components is then 

c = - 1 . (40) 

n na ' 
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For these computations the grid length, 21, was set equal to the wavelength, 
\, of the disturbance; so for the disturbance (primary) mode, n = 1, and 
in the linear range we can identify c^ with&p^ from the linear theory. 
Figure 10 shows that c^ is, indeed, a constant and the slope of this 
phase versus time plot shows that c^ = .363 which is in excellent 
agreement with Thomas' value of (K (3 = .3653 and. with our value of 

(k o xl = 3593- 

Other short computations were carried out with, again, a = 1.0, 

R^ = 6667, e = .005, hut with the initial mode shape, f^(y,0), chosen to 

correspond to one of the stable modes computed from linear theory. These 

calculations were carried out by arbitrarily setting f^(y,0) equal to 

* 

the mode § 2 g, and § 22 which are shown in Figure 9* The energy 

decay rates for these three cases are shown in Figure 10. From the energy 
decay rates and from the computed phase velocities (not shown here) the 
numerical integration technique yields the following values for the 
appropriate eigenvalues 

P 17 =1.351+i .138, P 2 q = 1.425 + i .207 

P 22 = 1.456 + i .0764. 

The linear matrix solution described in the preceeding section gave the 
values 

P 1? = 1.363 + i .136; P 2 g = 1.403 + i .192; 
p 22 = 1.462 + i .0742. 

* Since equation 33a is symmetrical in y it is possible to divide the 
solutions for $ nm (y) into even and odd solutions . The boundary conditions 
are §mn(+l) = $im(+l) = °> with $nm(°) “ (°) = 0 for even solutions and 

$nm(0) = $rlm(0) = ® for °dd solutions. For our numerical model with N 
lateral stations there are N-2 unknown (^jOnm'sj hence, N-2 eigenvalues 
are computed for each n, a, R„. Of these N-2 eigenvalues, N-l correspond 
to even solutions and N-3 correspond to odd solutions. ^ 



A more critical test of the numerical integration is the ability to 
predict the correct mode shape, § nm (y), for a disturbance of known wave — 
length and at a fixed R^ . Suppose in equation 38 we choose a - 1.0 
but arbitrarily select the initial disturbance mode shape f^(y,0). The 
arbitrary function f^y^O) can be expressed in terms of the eigenfunctions, 
§^(y), of the Orr-Sommerfeld equation (33a) for the chosen a and R^ 
(Schensted, 1961) . i.e., 



^(y.o) 



z 

m=l 



m 



V y) 



(Ula) 



with 



and 



A m ' ^ f l (y > 0) *m (y) ^ 



*m (y) * (72 - 
ay 




r*~> 



lm 



(41b) 



where $^(y) the solution to the adjoint of equation 33a. However, 
the linear solution tells us that only one of the in equation 4la is 
unstable; all other modes will decay with time. Hence, in the linear 
range only one mode, ^^(y), grows and will eventually dominate all the 
other terms in the expansion (4la). Therefore, the mode shape 



f-jyjt) - $ xl (y) as t - « 

where (y) is as shown in Figure 3 for the case a = 1.0, R = 6667. 
The initial shape chosen for q(y,0) is shown in Figure 11. The 

A 

disturbance kinetic energy, E(t), is given in Figure 12 and the mode 
shapes, f^(y,t), are shown for various times during the calculations. 
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As expected, the disturbance energy at first rapidly decreases then 
grows according to linear theory while the computed mode shape is 
clearly tending towards the predicted shape of ^ 1 (y) . 

From the results discussed in this section and from calculations 
at other wavenumbers, a, we have found that the numerical integration 
scheme consistently gives results that are in excellent agreement with 
the linear theory. The agreement between the results of the numerical 
model and the predicted results for small amplitude disturbances is 
gratifying and gives some confidence in the non-linear calculations 
now to be discussed. 

b) Finite Amplitude Disturbances 



Before discussing the results of the non-linear aspects of this 



investigation it is appropriate to give a very brief outline of the 
features of the existing non-linear theories which predict equations 
of the form of equation 1. The following discussion is from the 
paper by Pekeris and Shkoller (1967) which in turn outlines Stuart's 
(i960) theory for finite amplitude disturbances. If the Fourier 
representations for the disturbance streamfunction, Y' , and vorticity, 
Q ’ ; (equations 14 and 15) are substituted into the equation of motion 
(5a) the following equation for the modes results: 




tT 



O p f I! Hi *1 

- (nor) g n ) + i n a |^(U - f Q )g n - (U -f Q ) fnj (42) 



ia H 



n 



n = 0,1,2, . . . 



where 




(15c) 
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and 



H 



n 



Z 

s=l 



{(n-s)(f s ' g n _ s - g s ' 

y. f 

+ (n+s)(f g - 
v ' v s & n+s 



f 

n-s 



) 





(42a) 



If the parameter, | , is small (implying proximity to the 

neutral curve) then Stuart showed that in an approximation in which 
terms of order are retained, only the f Q , f^, and f^ in (l4a) 

are significant. Now expand f^(y,t) in terms of the eigenfunctions, 
of the Orr -Sommer f eld equation as was done in equation (4la) , 
so that 



yy,t) = e o l a (t) t (43) 

m=l 

Stuart’s asymptotic analysis for small e Q shows that only the first 
term in (43) is important 5 i.e., 

^(yjt) « e o A^t) $ lx (y) + 0 (e o 3 ) (44) 

and that the terms f and f„ are of smaller order so that 

o 2 

f 0 (y,t) - « 0 2 |A 1 (t)| a o (y) + 0(« o U ) (U 5 ) 

f 2 (y.t) - s Q 2 A^ft) G 2 (y) 0(e o ) lt (46) 



The functions G and G„ are determined from Pekeris and Shkoller's 
o 2 

equations ( 22 ), ( 23 ), and (24). By substituting equations (44), (45), 

and (46) into (42) one can obtain an expression of the form (l) for 

2 

the amplitude parameter, |A-^(t) | , i.e., 
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( 1 ) 



2 4 

— * S. + k 2 Kl 

The important point from our viewpoint is that according to the above 
analysis 5 the disturbance energy should be contained primarily in the 
f-L(y ? t) mode and that the growth of this energy should behave according 
to equation (l). Assuming that the disturbance streamfunction can be 
adequately represented by f^ then we have from (44) 

y' (x,y,t) = e Q (A 1 (t) § u (y)e‘ 10,x + A 1 (t) ^(y )e +1 (47) 

The disturbance kinetic energy as defined by equation (12c) (since we 
are neglecting the contribution from f^ the "disturbance" and "turbulent" 
kinetic energies are the same) becomes 

E (t) = e Q 2 I (A^t) | 2 (48) 

where 

1 = t 0 { I $ ii I 2 + 2 l $ nl 2 } ( u 9) 

It was pointed out in part 2 of this paper that the constant, k^, in 
(l) is related to the amplification factor of linear theory (k^= -2j 
and if k^ < 0 then the growth of an unstable disturbance is limited 
so that 

K f - ( 50 ) 

Pekeris and Shkoller found that these conditions are satisfied for a 
and R g inside the lower portion of the neutral curve. For example, 
at R e = 6667 and a - . 90 their calculations give 
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K \ 2 -* 1/9 

Oar solutions to the linear eigenvalue problem at R = 6667 and 

e 

a = .875 gives p = -.00394 and I = 1.848. Hence, 

k 

E(t) - (1.848) (.00394)/9.0 = 8.08 x 10 -4 

Figure 13 shows the energy growth for modes E Q , E^, E 2 , and E^ for 
a run with the initial disturbance given by (38) with o' = .875, 

R g = 6667, e = .005, and f^(y,0) = $-u_(y) where ^ ^ in the unstable 
eigenfunction for a = .875* The disturbance is, indeed dominated by 
E^(t) and the mode shape, f^(y,t), normalized so that f^(0,t) = 1.0, 
never departed from the initial shape given by §-^(y). A longer 
calculation was made using the same initial conditions except in this 
case a = 1.0. Again the energy remained predominately in the n = 1 
mode and again the mode shape for f^(y,t) never departed from the initial 
§11 (y) shape indicating that (44) and (47) do indeed adequately represent 
the flow as far as the overall disturbance energy is concerned. 

Although this latter calculation required 34 . 5 hours of computer 
time, we can see that the energy growth is still within the linear 
range. This demonstrates that, as was pointed out in part 2, it is 
impractical using the present techniques to follow the disturbance 
through the linear and into the non-linear range of amplitudes . For 
this reason we have been unable to establish whether or not the . 
"supercritical" equilibrium states as predicted by (l) and (50) do, 
in fact, exist for our exact treatment. Our calculations, as shown 
in Figures 13 and l4, do indicate, however, that there is no significant 
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difference between cases for which a = .875 or a = 1.0 in equation (38). 

The calculation of Pekeris and Shkoller, however , show that "super- 
critical” equilibrium states are possible for o' = .875 "but not for 
a = 1 . 0 . All our calculations for both low and high amplitude disturbances 
have shown no strong dependence on a as far as the final result is 
concerned. 

Although this last run was entirely within the linear energy 
growth range we will see that there were very interesting non-linear 
features to the calculations. The appearance of energy in modes 
E Q , E^, and E^ as shown in Figure l 4 is, of course, due to the 
presence of the non-linear terms in equation ( 5 a). In identical 
calcinations for which the non-linear terms in (5a) were deliberately 
omitted the primary (n=l) mode grew according to linear theory while 
the energy in the other modes remained constant or decreased and the 
energy levels (due to numerical "noise”) were twelve to fifteen orders 
of magnitude lower than that of the primary mode. Figure l 4 shows 
that the energy in the other modes (n = 0,2,3) grows very rapidly at 
first and then seems to follow an exponential growth rate, i.e., 

2p t 

E n (t) ~ e n n = 0 , 1 , 2, 3 

I 

for large t. From the data shown in Figure l 4 the values for P n 
are = . 0114 , ^ = .0050, = .0080, and = .0117. 
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The most interesting feature is the shape of the modes, f^, f^, 

* 

with respect to y . These shapes are shown in Figure 15 where the 

modes are normalized so that f (CLt) = 1.0 for even modes and for 

n 

odd modes the maximum real and imaginary parts of f (y) are unity. 

The phase velocities for these modes averaged between t = 65 and 
t - 71 were c^ - .362 and c^ = . 363 . These phase velocities are the 
same as the phase velocity of the primary mode predictable from linear 
theory (i.e^J ^ = . 363 ). 

The striking feature is the strong resemblence between these modes 
and the modes calculated from the Orr-Sommerfeld equation. Indeed, it 
appears that the non-linear solution has spontaneously generated modes 
that can be predicted from the linear theory even though the linear 
solution indicates that these are stable modes and would decay with 
time. This means that if we represent each mode by an expression of 
the form of (43) i.e., 

f (y,t) = % A (t) $ (y) (51) 

n 9 m=l m nm 



*From equation (15c) we see that g n has the same symmetry or anti- 
symmetry as does f n ; i.e. if f n (y) is an even function of y then g n 
is also even. Then from equation (42) we can establish that if 
Y'(x,y,0) is given by ( 38 ) and fq(y,0) is even in y then f Q (y,t) is 
odd, f 2 (y,t) is odd, f 3 (y,t) is even, f4(y,t) is odd, etc. Likewise 
if fq(y,0) is odd then fo(y,t) is odd, f 2 (y,t) even, f 3 (y,t) odd, 
^(yjt) is even, etc. For example, let fi(y,0) be even. Then since 
f n (y,0) = 0 n t 1 



3g 2 

dt 



-i aH, 



t=0 






t=o 



h " s i 



r i] 



= odd function of y. 
At the next time step, t = 6t 



g 2 (y 5 5t) =^ s 2/bt) 6t 



= odd function of y 
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then one mode seems to dominate the expansion. It is interesting to 
notice that in each case the dominant mode is one for which there exist 
a step gradient near the wall. Most of the stable modes that we have 
computed from linear theory are similar to the ones shown in Figure 9 
which have a detailed structure near the center of the channel rather 
than near the wall. 

A series of higher energy runs with various initial disturbance 
amplitudes was made at R g = 6667 using as initial conditions a 
disturbance with wave number a = 1.0 and that mode shape given by Thomas 
(see Figure 3)* The total energy growth (or decay) is shown in Figure 16 
for several different initial energy levels . A further energy breakdown 
is given in Figure 17 for the cases e = (.05)(/~5. For this case the 
total energy remains almost constant while there is a continual 
exchange of energy between the mean flow and the turbulence. The 

A 

turbulent kinetic energy oscillates about a value of E(t) = .03 which 
is roughly ten times the equilibrium value predicted from equation 50 . 

From Figure l4a it seems that there is a tendency towards some common 
energy equilibrium state for the different initial conditions . A more 
detailed inspection of the flow structure shows that there are significant 
differences between the three flows. Starting with the high energy run 
we summarize the structure of the flow by presenting plots of the mean 
velocity profile, Reynolds stresses, and energy spectrum at various 
times throughout the run (Figures 18, 19, and 20). Also, in Figure 21 
the shape of the primary (n = l) mode at the end of the calculations 
is shown for each of the three energy levels . For the high energy run 
(although the time parameter is comparatively small) there are drastic 
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alterations of the mean velocity profile and primary mode shape and 

considerable energy transfer from the primary mode to the higher 

harmonics. For the two lower energy runs there is little alteration 

of the parabolic velocity profile and only small amounts of energy 

transfer to higher wave numbers . For the lowest initial energy 

level (e - .05 / ^ 2) the primary mode shape remains essentially unchanged 

throughout the entire calculation. More details of this low energy 

run are shown in Figure 22 where the growth of the total turbulent 

energy in the first three spectral components (equation 19b) is 

shown. The phase velocities (equation 40) were computed near the 

end of the run and are summarized below: 

c ± = .492, c 2 = .484, c 3 = .473 

The measured C n were independent of y and all are roughly 40$ higher 
than the wave velocity (#tf3-^) the primary unstable mode predicted 
by linear theory. This low energy run has some qualitative similarities 
to the predictions of the non-linear theory which was outlined at the 
beginning of this section. That is, there appears to be an asymptotic 
energy level, the primary mode shape is only slightly altered and the 
energy in the primary mode remains dominant. 

In order to investigate the influence of initial wavenumber on 
the flow two runs were made at R g = 6667 with a - .875 (near the 
lower branch of the neutral curve) and with a = 1.05 (near the upper 
branch). For both cases the initial mode shape was determined from 
the eigenvalue solution to the linear equations and the initial 
energy was chosen to be slightly lower than that of the run (e = .05/ .f' 5> 
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Figure 1 6) which seemed to remain in an energy equilibrium state. The 
energy variation for these cases is shown in Figure 23. There is no 
obvious qualitative difference between the three wavenumbers; the initial 
high energy level, however, apparently causes oscillations of the 
disturbance kinetic energy and makes any direct comparison with non- 
linear theory rather difficult. 

Another important question which arises in the consideration of 
non-linear aspects in the stability of plane Poiseuille flows is whether 
or not a flow which is stable to infinitesimal disturbances might be 
unstable to disturbances of some finite amplitude. Such instabilities 
are generally referred to as "subcritical" instabilities . To investigate 
this aspect of the problem several runs were made using as initial 
conditions a disturbance that is known to be stable in the linear 
range. The case chosen was for a = .78, = 6667 for which p = .249 

+ i .00124 for the least stable eigenfunction determined from the linear 
solution described in Section 5* Again the initial disturbance was 
given by equation 38 with the function f^(y,0) = ^^(y) determined 
from the linear solution. For a small amplitude disturbance (e = .005) 
the disturbance energy decreased with time according to linear theory 
while for a larger amplitude (e - .0635). The energy variation was 
as shewn in Figure 24. The familiar oscillation of the turbulent 
kinetic energy which is associated with an initially high disturbance 
energy is present in this case; while the total kinetic energy appears 
to remain approximately constant after an initial decrease. 

We now address ourselves to the question posed at the end of 
Section 2 concerning the character of a two dimensional "turbulence” 
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in channel flows. For all the runs discussed above the computational 
grid length (2L) was equal to the wavelength of the initial unstable 
disturbance (i.e. 2L = 2tt/cv). In any finite difference representation 
of the equations of motion there are only discrete wavelengths 
available to be excited by the non-linear interactions . Recall from 
equation 24a that the streamfunction is expressed as 



Y 



K,J 



M 

Z 

n-1 



n, J 




(24a) 



and since Y v T is real 
J 




(Xf 



M-n+2, J 



^ f n,J "d f m-n+2, J 

i.e.,(£ f T is even and 'sj f is odd over the interval n = l,2,...M/2, 
. . .M. Therefore, there are M/2 distinct Fourier modes "available" to 
the solution of the equations of motion. When the grid length equals 
the wavelength of the initial disturbance there is only one lower mode 
available to be excited; namely, the zero mode, f . It is known 

0 3 J 

that two-dimensionality constrains the transfer of energy between 
the modes of which the disturbance is composed. For an isotropic 
two-dimensional turbulence there can be no systematic transfer of 
energy from lower to higher wave numbers without a corresponding 
transfer to still lower wavenumbers (Lilly (1968), Kraichnan (1967)). 
Hence, one may suspect that our discrete representation in which 
there are no available modes at lower wave numbers could force the 
disturbance energy to remain in the primary initial mode. The first 



alternative to this dilemma would seem to he to increase the grid 
length, 2L, thereby allowing modes with both higher and lower wave 
numbers to fit into the available grid length. However, if the 
initial disturbance energy is contained in, say, wave number n then 
the non-linear interaction of this mode with itself will initially 
excite wave numbers 0 and 2n; later modes 0, 2n, 3n, 4n, etc. The 
modes that can be excited will be separated by a spacing of n in 
wave number space, resulting in a spectrum with ,T holes ,T or regions 
in which no energy can appear. Clearly, this is not a realistic 
case. The alternative to this situation is to initially perturb the 
flow with a disturbance in which at least two contiguous (in the 
discrete wave number space) modes are present. For example, if the 
initial disturbance contains modes n and n+1 then the non-linear 
interaction will, at first, excite modes 0, 1, 2n, 2n+l, 2n+2; these 
in turn will excite modes 0, 1 , n-1, 2n-l, 4n, etc. until all possible 
wave numbers have been excited. 

It seems reasonable to require that the wavelengths of these two 
initial modes should both fall within the unstable region in the 
(a-R e ) plane. If we set the grid length so that 2L=16tt (or a = l/8) 
then modes with wavenumbers 7a and Qa will both fall within the 
unstable region. Hence 
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A lengthy numerical integration of the vorticity equations using 
the initial conditions just described is summarized in Figures 25 
and 26 which show the turbulent kinetic energy time variation and 
the energy spectrum at various times during the run. For this run 
the turbulent kinetic energy initially grows and then appears to 

A 

vary randomly about a value of E(t) = .019. The energy spectrum is 
quite different from the other cases presented in that there is a 
gradual transfer of significant amounts of energy to other wavenumbers. 
Although it is not shown here, the mean velocity profile again differed 
only slightly from the parabolic distribution. Although, there 
are significant differences between this run and the others presented, 
all the examples that have been discussed are qualitatively similar 
in that most of the energy remains in modes whose wavelength 
approximately equals that of the primary (unstable) mode. 

7. Conclusion 

Our numerical technique for integrating the vorticity equation 
has been shown to be stable and predicts results that are consistent 
with linear solutions to the vorticity equation. 

The influence of the finite difference representation on the 
linear eigenvalue solution to the equations of motion was investigated 
and some interesting conclusions resulted. The size of the 6y 
increment significantly influences the stability of the equations 
and for grids with 6y > .05 (approximately) the equations are stable 
to a .1 1 disturbance. Because of the DuFort-Frankel representation 
of the diffusion term, the growth rate of an unstable disturbance is 
strongly influenced by the 6t increment. The growth rate increases 
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rapidly with 6t even though the time step increment may still be 
small enough to satisfy ordinary numerical stability requirements . 

The 5x increment has the smallest effect on the linear solution 
and seems to be significant only when 6x is approximately equal to 
the channel half -width. 

Because of the long computing time required we were unable to 
follow the growth of an unstable disturbance through the linear 
and into the non-linear range of disturbance amplitudes . This 
precluded any direct comparison of our results with those predicted 
by the non-linear theories of Stuart and Watson and Reynolds and Potter. 
We can, however, make some general comparisons. We find qualitative 
agreement in the sense that the unstable primary mode dominates the 
disturbance and there is little distortion of the mean flow. The 
calculations presented indicate the existence of a disturbance energy 
equilibrium state which seems to be common (for a given Reynolds 
number) to all unstable disturbances despite their initial energy 
or wavenumber. 'Indeed* we have found no outstanding differences 
in the behaviour of disturbances with varying initial wave numbers 
(at a constant R ) which is not in agreement with the non-linear 
theories. The distortion of the primary mode shape is significant 
when the initial energy is high enough to be well beyond the linear 
range but for lower initial energies the distortion of the mode 
shape is less important. 

The computed flows which have been discussed in this paper 
show little resemblance to actual turbulent shear flows. In a true 
three dimensional flow one would expect an initial unstable disturbance 
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with a length scale approximately 2rr times the channel half width to 
grow in amplitude and through non-linear interactions the energy of 
the disturbance would be transferrred to higher wavenumbers with length 
scales comparable to the channel half width. The corresponding 
situation for a two dimensional turbulent channel flow is not clear. 
However, for two dimensional homogeneous turbulence the energy 
would cascade to the lower wavenumbers while mean square vorticity 
is transferred to higher wavenumbers when the turbulent energy is 
continually fed into the flow in a narrow band of spectral components. 
Neither of these processes was observed for the present calculations. 

One can speculate that our calculated "quasi -equilibrium" states 
in which most of the energy remains in one or two spectral components 
with little alteration of the parabolic velocity profile represent 
two dimensional solutions to the equations of motion which are, in 
some sense, highly unstable flows and, for this reason, are not 
observed experimentally. 
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Figure 2. Finite Difference Mesh 
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Figure 3. Comparison of the Eigenfunction 
5 , for R =6667 * a== 1 . 0 . 
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Figure 4. Effect of Grid Parameter, 6y, on 
the Disturbance Growth Rate 
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Figure 5. Effect of Grid Parameter, 
on the Eigenvalue, 
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Figure 6. Effect of Grid Parameter, 6t, on the 
Disturbance Growth Rat? 
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Figure 7. Energy Growth for a Low Amplitude Run. 

R = 6667, a = 1.0 
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Figure 8. Phase, cf^, Variation for a Low Amplitude 

Run . R = 6667 , a = 1 . 0 
e 
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Figure 9 a. Eigenfunction, $^(even), for R g = 6667 




s • 



Figure 9b. Eigenfunction, § 2 g(odd), for R e - 6667, 




o 




'■ 1 



Figure 9 c. Eigenfunction, $ (odd), for R = 6667, 






Figure 10. Computed Energy Decay for (a) (b) 

(c) $ 22 ’ Dashed line represents decay rate 
predicted from linear theory 
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Figure 11. Primary Mode Shape, f - » for a Low Energy 
(Linear) Run at Times t : P and t-6S. 
Dashed Line Represents £> u 
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Figure 12. Energy Variation for the Mode, f^(y,0), Shown 

in Figure 11. Dashed Line Represents Predicted 
Growth Rate for 
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Figure 13. Energy Variation for the Modes E o» E l«*2 
and E v 1^56667, c<* .375 
ss 





Figure 14. Energy Variation for the Modes E 0 , E^ , E2 
and E-,.K* C 6667, r 1.0 

5 * 



J?f 2 (y,t) _ /?f 2 (y,t) 




Figure 15a. Mode Shape for f Q (odd) at t = 85 . a = 1.0, 
R e = 6667. 2 

Dashed Line Represents $26 
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Figure 15b. Mode Shape for f 0 (even) at t = 85. or = 1.0, 

R e = 6667. 3 

Dashed Line Represents §^(even) 
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RELATIVE KINETIC ENERGY ( Total Minus Laminar) 








Figure 16. Total Kinetic Energy (Mean plus Turbulent) 
for Three Different Initial Energy Levels . 

The Laminar Energy (.600) has been Subtracted 
for Convenience. 

Top Curve, e = .5 

Middle Curve e = .05</lf" ( ot = 1.0, R =6667 
Bottom Curve e = .05 /a/2 J 
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Figure 17. Energy Variation for the Run e = .05 J~5 
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Figure l8a. Mean Velocity Profiles for the Run e = .5, 
a = 1.0, R e = 6667 

Dashed Line is U(y) = 3/2(1 - y 2 ) 
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Figure 18b. Reynolds Stresses for the Run e = . 5 , a = 1.0 
R = 6667 

e 
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Figure 19a. Mean Velocity Profiles for the 

Run e = .05\/5j ot = 1.0, R =6667 

^ p 

Dashed Line is U(y)=3/2(l-y ) 



GU 





A A 

U V 



65 



En (t) 




Figure 19c. Energy Spectrum for the Run e = .05*/ r 5~, 
a = 1.0, R = 666 7 at t = 59 
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Figure 20 a. Mean Velocity Profiles for the Run 

s = .05/1/27 a - 1 . 0 , R g = 6667 
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Figure 20b, Reynold* s Stresses for the Run e = .05//2’ 

a = 1.0. R = 6667 
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Figure 20c. Energy Spectrum for the Run e 

a = 1.0. R = 6667 69 
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Figure 21a. Primary Mode Shape, f , for the Run e = .5 
at t = 11 
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Figure 21b. Primary Mode Shape, f , for the Run 

e = . 05 JT at t = 59 
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Figure 21c. Primary Mode Shape, for the 
e = .05A/2 - at t - 112 





Figure 22. Energy Variation in E^, and E^ for the 
Run e = .05A/2~ 
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Figure 23 . Energy Variation at R 0 = 6667 for Three 
Different Wavenumbers . 

- 0 - Turbulent Kinetic Energy 
-x- Total Kinetic Energy 
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R e = 6667, e = .0635. Dashed 

line is the total kinetic energy 
while solid line is turbulent 
kinetic energy. 
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Figure 25. Turbulent Kinetic Energy for the Run with Modes 
f and f o Initially Present R = 6667, e = .059 
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Figure 26. Energy Spectrum for the Run with 
Modes f and fg Initially Present 

R = 6667, e = .059 
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